IMPC Mouse data - Variance in sex differences

# loads the raw data, setting some default types for various columns

load_raw <- function(filename) {
  read_csv(filename, 
                              col_types = cols(
                                .default = col_character(),
                                project_id = col_character(),
                                id = col_character(),
                                parameter_id = col_character(),
                                age_in_days = col_integer(),
                                date_of_experiment = col_datetime(format = ""),
                                weight = col_double(),
                                phenotyping_center_id = col_character(),
                                production_center_id = col_character(),
                                weight_date = col_datetime(format = ""),
                                date_of_birth = col_datetime(format = ""),
                                procedure_id = col_character(),
                                pipeline_id = col_character(),
                                biological_sample_id = col_character(),
                                biological_model_id = col_character(),
                                weight_days_old = col_integer(),
                                datasource_id = col_character(),
                                experiment_id = col_character(),
                                data_point = col_double(),
                                age_in_weeks = col_integer(),
                                `_version_` = col_character()
                              )
  )
}

# Apply some standard cleaning to the data
clean_raw_data <- function(mydata) {
  mydata %>% 
    
    # Fileter to IMPC source (recommened by Jeremey in email to Susi on 20 Aug 2018)
    filter(datasource_name == 'IMPC') %>%
    
    # standardise trait names
    mutate(parameter_name = tolower(parameter_name) ) %>%
    
    # remove extreme ages
    filter(age_in_days > 0 & age_in_days < 500) %>% 

    # remove NAs 
    filter(!is.na(data_point)) %>%
  
    # subset to reasonable set of variables
    # date_of_experiment: Jeremy suggested using as an indicator of batch-level effects
    select(production_center, strain_name, strain_accession_id, biological_sample_id, pipeline_stable_id, procedure_group, procedure_name, sex, date_of_experiment, age_in_days, weight, parameter_name, data_point) %>% 
    arrange(production_center, biological_sample_id, age_in_days)
}

# this is a necessary step for the loop across all traits 
data_subset_parameterid_individual_by_age <- function(mydata, parameter, age_min, age_center) {
  tmp <- mydata %>%
    filter(age_in_days >= age_min,
           id == parameter) %>%
    # take results for single individual closest to age_center
    mutate(age_diff = abs(age_center - age_in_days)) %>%
    group_by(biological_sample_id) %>%
    filter(age_diff == min(age_diff)) %>%
    select(-age_diff)
  
  # still some individuals with multiple records (because same individual appears under different procedures, so filter to one record)
  j <- match(unique(tmp$biological_sample_id), tmp$biological_sample_id)
  tmp[j, ] 
}

calculate_population_stats <- function(mydata, min_individuals = 5) {
  mydata %>% 
    group_by(population, strain_name, production_center, sex) %>% 
    summarise(
      trait = parameter_name[1],
      x_bar = mean(data_point),
      x_sd = sd(data_point),
      n_ind = n()
    ) %>% 
    ungroup() %>%
    filter(n_ind > min_individuals) %>% 
    # Check both sexes present & filter those missing
    group_by(population) %>% 
    mutate(
      n_sex = n_distinct(sex)
    ) %>% 
    ungroup() %>%
    filter(n_sex ==2) %>% 
    select(-n_sex) %>%
    arrange(production_center, strain_name, population, sex)
}

create_meta_analysis_effect_sizes <- function(mydata) {
  i <- seq(1, nrow(mydata), by = 2)
  input <- data.frame(
    n1i = mydata$n_ind[i],
    n2i = mydata$n_ind[i + 1],
    x1i = mydata$x_bar[i],
    x2i = mydata$x_bar[i + 1],
    sd1i = mydata$x_sd[i],
    sd2i = mydata$x_sd[i + 1]
  )
  
  mydata[i,] %>% 
  select(strain_name, production_center, trait) %>%
    mutate(
      effect_size_CVR = Calc.lnCVR(CMean = input$x1i, CSD = input$sd1i, CN = input$n1i, EMean = input$x2i, ESD = input$sd2i, EN = input$n2i),
      sample_variance_CVR = Calc.var.lnCVR(CMean = input$x1i, CSD = input$sd1i, CN = input$n1i, EMean = input$x2i, ESD = input$sd2i, EN = input$n2i),
      effect_size_VR = Calc.lnVR(CSD = input$sd1i, CN = input$n1i, ESD = input$sd2i, EN = input$n2i),
      sample_variance_VR = Calc.var.lnVR(CN = input$n1i, EN = input$n2i),
      effect_size_RR = Calc.lnRR(CMean = input$x1i, CSD = input$sd1i, CN = input$n1i, EMean = input$x2i, ESD = input$sd2i, EN = input$n2i),
      sample_variance_RR = Calc.var.lnRR(CMean = input$x1i, CSD = input$sd1i, CN = input$n1i, EMean = input$x2i, ESD = input$sd2i, EN = input$n2i),
      err = as.factor(seq_len(n()))
    )

}

Calc.lnCVR <- function(CMean, CSD, CN, EMean, ESD, EN){
  log(ESD) - log(EMean) + 1 / (2*(EN - 1)) - (log(CSD) - log(CMean) + 1 / (2*(CN - 1)))
}

Calc.var.lnCVR <- function(CMean, CSD, CN, EMean, ESD, EN, Equal.E.C.Corr=T) {
  if(Equal.E.C.Corr==T){
    mvcorr <- 0 #cor.test(log(c(CMean, EMean)), log(c(CSD, ESD)))$estimate   old, slightly incorrect
    S2 <- CSD^2 / (CN * (CMean^2)) + 1 / (2 * (CN - 1)) - 2 * mvcorr * sqrt((CSD^2 / (CN * (CMean^2))) * (1 / (2 * (CN - 1)))) + ESD^2 / (EN * (EMean^2)) + 1 / (2 * (EN - 1)) - 2 * mvcorr * sqrt((ESD^2 / (EN * (EMean^2))) * (1 / (2 * (EN - 1))))
  }
  else{
    Cmvcorr <- cor.test(log(CMean), log(CSD))$estimate
    Emvcorr <- cor.test(log(EMean), (ESD))$estimate
    S2 <- CSD^2 / (CN * (CMean^2)) + 1 / (2 * (CN - 1)) - 2 * Cmvcorr * sqrt((CSD^2 / (CN * (CMean^2))) * (1 / (2 * (CN - 1)))) + ESD^2 / (EN * (EMean^2)) + 1 / (2 * (EN - 1)) - 2 * Emvcorr * sqrt((ESD^2 / (EN * (EMean^2))) * (1 / (2 * (EN - 1))))     
  }
  S2
}

Calc.lnVR <- function(CSD, CN, ESD, EN){
  log(ESD) - log(CSD) + 1 / (2*(EN - 1)) - 1 / (2*(CN - 1))
}

Calc.var.lnVR <- function( CN,  EN) {
  1 / (2*(EN - 1)) + 1 / (2*(CN - 1))
}

Calc.lnRR <- function(CMean, CSD, CN, EMean, ESD, EN) {
  log(EMean) - log(CMean)
}

Calc.var.lnRR <- function(CMean, CSD, CN, EMean, ESD, EN) {
  CSD^2/(CN * CMean^2) + ESD^2/(EN * EMean^2)
}

## Load raw data - save cleaned dataset as RDS for reuse
#data_raw <- load_raw("data/dr7.0_all_control_data.csv") %>%
#    clean_raw_data()
#dir.create("export", F, F)
#saveRDS(data_raw, "export/data_clean.rds")
getwd()

data1 <- readRDS("../export/data_clean.rds") 
#data1 

data3 <- data2 %>%
mutate(population = sprintf("%s-%s", production_center, strain_name))

 group <- read.csv("../export/ParameterGrouping.csv") 
 data <- data3
 data$parameterGroup <- group$parameter[match(data$parameter_name, group$parameter_name)] 

#head(data)

names(data)[16] <- "parameter_group"    

data <- transform(data, id = match(parameter_name, unique(parameter_name)))
n1 <- length(unique(data$parameter_name)) #232
n2 <- length(unique(data$parameter_group)) #161
n3 <- length(unique(data$procedure_name)) # 26

n <- length(unique(data$id))
#n  # just to check that the number of traits is 232

results.alltraits.grouping <- as.data.frame(cbind(c(1:n), matrix(rep(0, n*14), ncol = 14))) #number of individual results per trait = 10
names(results.alltraits.grouping) <- c("id", "lnCVR", "lnCVR_lower", "lnCVR_upper", "lnCVR_se", "lnVR", "lnVR_lower", "lnVR_upper", "lnVR_se", "lnRR", "lnRR_lower", "lnRR_upper" ,"lnRR_se" , "sampleSize", "trait")

meta1 <- meta_clean 
# length(unique(meta1$procedure)) #18
# length(unique(meta1$GroupingTerm)) #9 
# length(unique(meta1$parameter_group)) # 148 levels. To be used as grouping factor for meta-meta analysis / collapsing down based on things that are classified identically in "parameter_group" but have different "parameter_name" 
length(unique(meta1$parameter_name)) #218

kable(cbind (meta1 %>%  count(parameter_group) )) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "200px")

meta1b <-
  meta1 %>%
  group_by(parameter_group) %>% 
  summarize(par_group_size = length(unique(parameter_name, na.rm = TRUE)))
#this gives a summary of number of parameter names in each parameter group, now it neeeds to get merged it back together


meta1$par_group_size <- meta1b$par_group_size[match(meta1$parameter_group, meta1b$parameter_group)]

# Create subsets with > 1 count (par_group_size > 1) 

meta1_sub<-subset(meta1,par_group_size >1) # 90 observations   
meta1_sub$sampleSize <- as.numeric(meta1_sub$sampleSize)

count_fun <- function(mod_sub)
  return(c(mod_sub$reg_table$b.r, mod_sub$reg_table$CI.L, mod_sub$reg_table$CI.U, mod_sub$reg_table$SE))   #estimate, lower ci, upper ci, SE

robusub_RR <- model_count %>% 
  transmute(parameter_group, estimatelnRR = map(model_lnRR, count_fun)) %>% 
  mutate(r = map(estimatelnRR, ~ data.frame(t(.)))) %>%
  unnest(r) %>%
  select(-estimatelnRR) %>%
  purrr::set_names(c("parameter_group","lnRR","lnRR_lower","lnRR_upper","lnRR_se"))

robusub_CVR <- model_count %>% 
  transmute(parameter_group, estimatelnCVR = map(model_lnCVR, count_fun)) %>% 
  mutate(r = map(estimatelnCVR, ~ data.frame(t(.)))) %>%
  unnest(r) %>%
  select(-estimatelnCVR) %>%
  purrr::set_names(c("parameter_group","lnCVR","lnCVR_lower","lnCVR_upper","lnCVR_se"))

robusub_VR <- model_count %>% 
  transmute(parameter_group, estimatelnVR = map(model_lnVR, count_fun)) %>% 
  mutate(r = map(estimatelnVR, ~ data.frame(t(.)))) %>%
  unnest(r) %>%
  select(-estimatelnVR) %>%
  purrr::set_names(c("parameter_group","lnVR","lnVR_lower","lnVR_upper","lnVR_se"))

robu_all <- full_join(robusub_CVR, robusub_VR) %>% full_join(., robusub_RR)
kable(cbind(robu_all, robu_all)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "200px")

meta_all <- meta1 %>% filter(par_group_size == 1) %>% as_tibble
#str(meta_all)
#str(robu_all)
#which(is.na(match(names(meta_all),names(robu_all))))  # check

# Step1 
combinedmeta <- bind_rows(robu_all, meta_all)
#glimpse(combinedmeta)

# Steps 2&3
metacombo <- combinedmeta
metacombo$counts <- meta1$par_group_size[match( metacombo$parameter_group, meta1$parameter_group)]
metacombo$procedure2 <-meta1$procedure[match( metacombo$parameter_group, meta1$parameter_group)]
metacombo$GroupingTerm2 <-meta1$GroupingTerm[match( metacombo$parameter_group, meta1$parameter_group)]

kable(cbind (metacombo, metacombo)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "200px")